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Abstract 

A one-dimensional model evolution equation is used to describe the nonlinear dy- 
namics that can lead to the breakup of a cylindrical thread of Newtonian fluid when 
capillary forces drive the motion. The model is derived from the Stokes equations by 
use of rational asymptotic expansions and under a slender jet approximation. The 
equations are solved numerically and the jet radius is found to vanish after a finite 
time yielding breakup. The slender jet approximation is valid throughout the evolu- 
tion leading to pinching. The model admits self-similar pinching solutions which yield 
symmetric shapes at breakup. These solutions are showm to be the ones selected by 
the initial boundary value problem, for general initial conditions. Further more, the 
terminal state of the model equation is shown to be identical to that predicted by a 
theory which looks for singular pinching solutions directly from the Stokes equations 
without invoking the slender jet approximation throughout the evolution. It is shown 
quantitatively, therefore, that the one-dimensional model gives a consistent terminal 
state with the jet shape being locally symmetric at breakup. The asymptotic expan- 
sion scheme is also extended to include unsteady and inertial forces in the momentum 
equations to derive an evolution system modelling the breakup of Navier-Stokes jets. 
The model is employed in extensive simulations to compute breakup times for differ- 
ent initial conditions; satellite drop formation is also supported by the model and the 
dependence of satellite drop volumes on initial conditions is studied. 
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1 Introduction 


It is well known that, a perfectly cylindrical jet of liquid which supports surface tension at 
its free interface is an unstable stationary solution of the equations of motion (both viscous 
and inviscid) - any uniform axial velocity can be removed by a Galilean transformation, so 
we concentrate on stationary unperturbed jets. The instability is driven by capillary forces: 
If the radius of the jet decreases locally at some axial position, capillary forces induce a 
local increase in pressure (the pressure outside the jet is constant here) and conversely a 
local increase in radius causes a local decrease in pressure just below the interface; fluid will 
flow from regions of interfacial depression to regions of interfacial expansion and by mass 
conservation the depression/expansion will decrease/increase. Linear stability analyses based 
on normal modes have been carried out by Rayleigh for both inviscid (see [1]) and viscous 
(see [2] ) jets with the effect of the surrounding fluid neglected. The effect of a viscous 
surrounding phase has been included in the linear stability analysis of Tomotika [3]. The 
viscous dispersion relation given in [2] is valid in the limit of a highly viscous fluid; the 
general dispersion relation for arbitrary viscosities is given in [3] as well as Chandrasekhar 
([4], p. 541). The present work considers Stokes flows and so Rayleigh’s dispersion is 
useful, then. Following Tomotika, normal mode disturbances are considered proportional to 
exp(?(n/ + £■*)), where t is time, z is axial distance, k is the wavenumber of the disturbances 
and in their growth rate: Solution of the linear eigenvalue problem in the case of Stokes jets 
with the surrounding phase neglected, yields the following dimensional growth rates (see [3], 

[ 4 ]): 


(7 k 2 R 2 - 1 

2Rji k 2 R 2 + 1 - k 2 R 2 I 2 (kR)/I 2 (kR)' 


( 1 ) 


where u is the surface tension coefficient (assumed to be constant), ft is the fluid viscosity 
and R is the unperturbed jet radius. According to (1), the most unstable wave has k = 0 
and the growth rate is then given by 


in o 


<7 

6 Rfi 


(2) 


This long-wave instability result is due to the fact that time derivatives are dropped in 
the momentum equations. Its inclusion (or inclusion of a surrounding phase) provides a 
dispersion relation with a non- zero maximally growing wavenumber - see Tomotika’s results. 
Our interest in this article is to generate analytical nonlinear structures which coincide with 
boundary integral numerical solutions, for example, in the long wave regime and in particular 
near pinching. The theory presented in this article, then, is a nonlinear long wave one and 
so on linearization of the obtained evolution equations, the growth rates (2) emerge (see 
below). It is worth noting that since the maximum growth rate is for infinitely long waves 
(as predicted theoretically by linear theory at least), a long wave nonlinear ansatz is justified. 
For inviscid and viscous (but not Stokes) jets the maximum growth rate occurs at a finite 
wavenumber which lies between 0 and 1 (in non-dimensional terms) and long wave theories, 
then, may not be consistent with the linear results; they can, however, describe phenomena 
such as jet pinching, which are beyond the scope of linear theory (see later). 

Experimental observations of jet breakup phenomena have been carried out by Chaud- 
hary and Maxworthy [5], [6], Donelly and Glaberson [7], Goedde and Yuen [8] and Peregrine, 
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Shoker and Symon [9]; the latter investigation considers breakup with gravity being impor- 
tant also. More recently Tjahjadi, Stone and Ottino [10] have carried out an experimental 
study with highly viscous fluids using a Couette device; numerical boundary integral simu- 
lations based on the Stokes equations are also described and compared to the experiments. 
This work is an extension of previous experiments and computations by Stone and Leal [11] 
who consider the breakup of an initially extended drop. Simulations based on the Navier- 
Stokes equations and different Reynolds numbers (covering inviscid to viscous flow’s) have 
recently been applied by Richards, Lenhoff and Berris [12] to describe the motion of a vis- 
cous jet injected into another viscous liquid. Comparisons of the simulated breakup with 
experiments is made with good overall agreement. All the simulations described are carried 
out for axisymmetric jets; this assumption is consistent with linear theory as well as many 
experimental situations especially if the Reynolds number is not too large. For a review of 
drop formation in circular jets, see Bogy [13]. 

Our interest is in the description of the breakup from the viewpoint that it is a singu- 
larity of the equations of motion with the jet radius vanishing after a finite time at some 
axial location. The theory developed is a local one and the singular time as well as the 
axial position where breakup first occurs, depend on initial conditions. The main reason 
why such local structures are desirable analytically is that they provide a rational way of 
continuing the solutions beyond breakup after the topology changes. In fact simple mass and 
moment um balances can provide regular initial conditions for the dynamics beyond breakup 
by assuming that a spherical blob of fluid is attached to the end of the breakaway jet (see 
Ting and Keller [14] for inviscid flow's and Papageorgiou [15] for Stokes and Navier-Stokes jet 
flows). Our approach is based on working with a simplified set of evolution equations which 
arise from the governing equations after a slender jet ansatz is adopted. One-dimensional 
slender jet models have recently been used by many authors in the modeling of viscous liquid 
jets (see for example Renardy [16], Eggers and Dupont [17], Eggers [18], [19], and Garcia 
and Castellanos [20]). Such approaches are unsteady extensions of the steady fiber extru- 
sion problem considered by previous investigators (see for example Schultz and Davis [21] 
and references therein, as well as the review article by Denn [22]). In the present w’ork w T e 
mostly consider breakup governed by the Stokes equations (extensions which include inertail 
and unsteady terms in the momentum equations are also derived); this is a problem with 
applications in microgravity flows, for instance, where fluid viscosities are high and typical 
Reynolds numbers are small. A system of one-dimensional partial differential equations is 
derived by assuming that the ratio between the maximum jet radius to the wavelength of in- 
terfacial deflections is an asymptotically small parameter w'hich can be used in an asymptotic 
expansion to capture the leading order evolution. (Note that, higher order corrections can be 
calculated routinely within this framew'ork). These equations have been derived previously 
in [16] using a Lagrangian coordinate system. Renardy also proves a theorem which states 
that in the Newtonian case these equations have solutions with the radius vanishing after 
a finite time (he also show's that this is not the case for several viscoelastic models). The 
equations support similarity solutions corresponding to pinching and w T e show' that these 
are identical to pinching solutions obtained directly from the full Stokes equations (see also 
[15] where the Navier-Stokes system is also studied). Analysis of the similarity equations 
fixes universal scaling laws for the breakup and we confirm these scaling exponents by direct 
numerical solution of the model equations for different, initial conditions. 
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The article is organized as follows. In Section 2 the equations together with interfacial 
boundary conditions are given along with the non-dimensionalization based on capillary 
scales. Section 3 derives the slender jet model comprised of a system of coupled nonlinear 
evolution equations. A generalization of the expansions to include inertial and unsteady 
momentum effects is given in an Appendix. In Section 4 we construct self-similar solutions 
valid as the jet pinches after a finite time. This is done for (i) the evolution equations derived 
in Section 2, and, (ii) by looking for singular solutions directly from the full Stokes equations 
given in Section 2. The similarity equations are the same and they are solved to obtain the 
solutions in closed implicit form. The scaling exponents are also determined and are shown to 
be universal. Section 5 is devoted to the numerical solution of the evolution equations. It is 
shown that the analytical self-similar solutions of Section 4 are the ones obtained at pinching. 
This is done by comparison of scaling exponents and functions provided by the simulation 
with the analytical predictions. Both symmetric and non-symmetric initial conditions (for 
the jet shape) are used and the analytical results are verified in both cases. 


2 Governing equations 


Consider the evolution of a viscous cylindrical column of fluid of viscosity p. Initially the 
jet is a perfect infinite cylinder of radius R and zero velocity (a constant axial velocity 
can be removed by a suitable Galilean transformation); this is an exact solution of the 
equations of motion and boundary conditions and is the flow used in linear stability theories 
in the calculation of the initial stages of the instability. Even though linear theory cannot 
provide a quantitative description of the flow at breakup, it does give useful insight into 
the competing physical mechanisms acting. It is found, then, that capillary forces drive the 
instability which leads to pinching (see Introduction). The pertinent scales are: lengths scale 
with /?; velocities scale with \ pressure scales with time scales wfith l ~. The equations 
are made dimensionless by introduction of the variables 

( r,z ) = R(r, z), ( u,w ) = -{u,w), p = ^ p , t - —t. (3) 

// rt (7 


Substitution of (3) into the Stokes equations and interfacial boundary conditions and drop- 
ping of the bars, yields the following non-dimensional system: 


Au -u = p T , 

Aw = p z , 

- {ru) + 10 z = 0, 

r 


(4) 

(5) 

( 6 ) 
(7) 


-t q 2 * t * # * 

where A = 4^ + + f: ?• The interfacial conditions of tangential stress balance, normal 

stress balance and kinematic condition, on r = S{t,z) are 


( u z + te r )(l - s 2 z ) + 2u t S z - 2w z S z = 0, 
p — 2u r — (— p + 2 w z )S 2 + 2 (u z + w r )S z — — ^S zz — — (1 + S 2 2 )) (1 + S^) ^ 2 , 

u = S t + wS z . 
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( 8 ) 

(9) 

(10) 



Iii addition to the interfacial conditions (S)-(10) we impose regularity of flow quantities on 
the jet axis r = 0. 

We are interested in solutions of the system (4)-(10) in the strongly nonlinear regime and 
in particular when pinching occurs and the jet radius tends to zero at some axial position. 
Direct numerical solutions have been given in [10] and [11] by use of boundary integral tech- 
niques. Here we propose a quantitative description of the pinching by use of an asymptotic 
theory. 

3 The evolution equations governing pinching 

The evolution equations arise from a nonlinear long wave theory. If we assume that a 
typical dimensional length scale in the axial direction is D , then the ratio e = ■§ is taken 
to be a small quantity and we seek a leading order solution of an asymptotic expansion 
in e. Such approximations have been used extensively in the description of nonlinear long 
wave interfacial flows (see Introduction). At first sight the theory seems arbitrary since the 
parameter D does not have a distinct physical meaning. The usefulness of the evolution 
equations given below along with their solutions lies in the fact that as a pinch forms the 
radial length scale is at most of the order of the axial length scale (i.e. the axial extent of 
the pinch region); if the radial length scale is asymptotically smaller than the axial one, a 
long wave approximation is valid to leading order and so the model equations describe local 
pinching solutions of the full equations even though the transient motion may be inaccurate. 
A theory which constructs pinching solutions of the Stokes equations directly, has been given 
in ([15]). It is shown later that the two theories are in complete agreement thus lending weight 
to the relatively simple model equation approach. 

The long wave ansatz is easily applied to (4)-(10) by introducing the transformation 
Jh — > ejh; in addition the appropriate expansions proceed in powers of e 2 (this can be seen 
from the streamfunction formulation of the Stokes equations in cylindrical coordinates which 
yields the biharmonic equation for the streamfunction implying that the velocities expand 
in powers of e 2 ): 


u(t, r, z ) — Uq T e 2 Ui -T . . . , 

(ID 

?(t, r, z) = -j-u’o + at’i + . . . , 

(12) 

p(t, r, z) = po + e 2 p\ T . . . , 

(13) 

S(t, z) = 5q + e 2 >?i + . . . . 

(14) 


Substitution of (1 1 )-( 14) into the governing equations and boundary conditions gives a se- 
quence of problems. The leading order velocity components are found from (5) and (6) and 
are 

IV 0 = w 0 {t,z), u 0 = (15) 

Substitution of this expression for u 0 into (4) gives p 0r = 0 and so p 0 = p 0 (t,z). This readily 
yields the correction uq from (5) as well as the leading order pressure distribution throughout 
the jet from the normal stress balance condition (9); the continuity equation then gives Ui in 
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terms of iv 0 and the process can be continued to higher order. The results which are needed 
for the evolution equation are: 

uq(f,r,2) = -\r 2 (w 0 zz -poz) + A(t,z), po(t,z) = — w 0z , (16) 

4 oo 

where A[t,z) is to be determined but it does not enter the leading order dynamics. With 
these solutions available, the evolution equations are obtainable from the tangential stress 
balance condition (8) and the kinematic condition (10). To leading order the tangential 
stress balance gives Wo r = 0 which is already satisfied by (15). Next at order e and at order 
one in (8) and (10) respectively, we find 

u>0z + w lr + 2tZ0 r iSos — 2-WozSoz = 0, r = So(t , ^), (1 0 

u o — Sot + w qSq z . (18) 

The desired evolution equations are obtained from (17) and (18) by elimination of u 0 and uq 
in favor of u'o alone to yield an evolution system for So and the leading order axial velocity 
uv This system is: 


Sot + woSq z + -wozSo — 0, (19) 

-Swoz 4- 3Soz w oz + = (20) 

Z Zoo 

Asymptotic and numerical solutions of (19) and (20) are given in later Sections but some 
comments on the physical origin of these equations are useful. 

An integrated form of the evolution equations has been derived previously by M. Renardy 
[16] by use of physical arguments. The first equation describes conservation of mass for 
slender jets which is most easily seen by multiplication of (19) by So to yield the conservation 
form 


(S|)< + (“-’oSo), = o. (21) 

The second equation was written down by Renardy by considering the force acting on a 
cross-section of the slender jet. The Stokes flow has no inertia and so this force must be 
constant along the jet; equation (20) is the ^-derivative of this force. This is most easily 
seen by multiplication of (20) by So and writing it in the form 

(3S 0 W + S 0 )z = 0. (22) 

These ideas have also been used in [15] in following jet evolution just beyond pinching by 
description of overall mass and momentum balance equations (see also [14] for an application 
to inviscid flows). 

We note also that a simple modification of (11 )-( 14) that allows for non zero Reynolds 
numbers and introduces inertial and unsteadiness into the momentum equations leads to 
a set of evolution equations which can be used to model pinching of jets governed by the 
Navier-Stokes equations. This analysis is included in Appendix A. 
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4 Finite-time singularities. Self-similar solutions. 

The evolution equations (19) and (20) have been proven by Renardy [16] to possess singular- 
ities with the jet radius vanishing after a finite time. A Langrangian formulation was used to 
prove the theorem and in particular it is established that arbitrarily small initial conditions 
can lead to breakup. Our interest is in the related problem of quantifying such terminal 
states and in particular in establishing any type of dominant solutions at breakup. To this 
end we employ numerical computations of the initial value problem and a local analysis of 
breakup by construction of self-similar solutions. In what follows we address the description 
of local structures at. breakup. This is done for (i) the long wave evolution system derived 
above, and, (ii) directly from the full Stokes equations. The latter analysis has been given in 
[15] where unique scaling laws were established by solution of a nonlinear eigenvalue problem 
(see later). In what follows we show that the similarity solutions of the model equations are 
identical to those found by a direct analysis of the Stokes system and we test the analytical 
self-similar structures and in particular the unique scaling exponents, with numerical solu- 
tions obtained by solving the initial value problem (19) and (20). This is done in a later 
Section with excellent agreement. 


4.1 Pinching solutions of the model equations 

We work with (19) and (20) and assume that the jet radius vanishes after a finite time, t s say. 
The objective is to describe the solutions near this time, so that r = t s — t with 0 < r << 1, 
and near the axial position where the jet breaks (without loss of generality this position 
is taken to be the origin). A balance of terms in (19), (20) gives the order-of-magnitude 
estimates 

So WqSo SqWo 

~ ~ ‘^ 0 * 

T Z z 

from which it follows that So ~ t while tco ~ Assuming that the phenomenon is a 
focussing one we can introduce a positive parameter ft which controls the extent of the 
similarity region by £ ~ t 13 . Formally, then, we look for pinching solutions of (19) and (20) 
in the form 

S 0 (M) = T .f{0, u-o {t,z) = T 3 ~ l g{(), £ = ^ (r = t 5 -t). (23) 


Note that the forms (23) above are an exact self-similar transformation since all terms in 
the governing equations are in balance. The situation is slightly different for the full Stokes 
equations (see later). Partial derivatives transform according to 


d_ _d_ fcd_ d_ ]_d_ 

dt dr ^ t <9£ ’ dz t& d £ 


(24) 


Substitution of (23) and (24) into (19) and (20) gives the following equations for the scaling 
functions: 


(s + K)f + (\s'-i)f = 0’ 

|(3/V + /)=0. 
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(25) 

(26) 



Equation (26) can be integrated once to give 


+ p (27) 

where k is a constant of integration. Equations (25), (27) need to be solved and the constants 
/3 and k determined. In fact the value of k can be expressed in terms of f as explained next. 
The force balance equation (22) can be integrated in z to yield (44) below and the function 
A(/) given by (45). £.From the ansatz (23) and (44) it is clear that A ~ r as r — + 0 and in 
fact A = 3 kr -f . . . then with k as used in (27) above. The following expression for k then 
follows, 

, _ 1 

from either of two equivalent ways: (i) by integration of (27) over the range of £ and using 
the fact that axial velocities are zero far from the pinch (this is shown asymptotically later), 
(ii) by introduction of the ansatz (23) into the expression (45) for A. This identification of 
k is essential in the determination of a unique value of (3 described later. Before doing this 
we present the construction of singular solutions of the full Stokes equations which yield an 
identical result . 

4.2 Pinching solutions of the Stokes equations. 

In this section we summarize results described in more detail in [15]. The idea is to construct 
solutions to (4)- (10) with the jet radius going to zero after a finite time. The system is not 
one dimensional in space as for the model equations and a similarity variable in the radial 
direction is also needed. The following transformations are appropriate 

r = r°y , z = S = t"/((), (29) 

w = T’’W(t,p,(), u = t' +c ’~ 0 U(I p = T~°P, (30) 

where u follows from the continuity equation and a > (3 in keeping with slender geometries 
at breakup. The constants a,/?, 7 are to be determined along with the scaling functions (the 
use of the same symbols as before should not be confusing). Unlike the model equations, 
the transformations (30) do not retain all terms in the Stokes equations and an asymptotic 
expansion in powers of r 2a ~ 213 is appropriate (this comes from the biharmonic operator for the 
streamfunction in much the same way as the e 2 expansion was established for the derivation 
of the model system). These expansions are 

«■’ = r* (<?(£) + T 2a ~ 2/) W 1 + ...), « = r^ a ~ 0 {-\yg 0 ( + r 2a ~ 20 U\ + ...), 

which on substitution into the z-momentum equation (5) and balance of leading order terms 
yields 

7 = j 3- Q , ir„ = \y(P< - «’«{). (31) 


(28) 
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This in turn leads to lh and the result that P y = 0. Substitution of leading order solutions 
into the tangential stress balance equation (8) and the normal stress balance (9) give at 
0{t~ P) and 0{r~ a ) respectively 




on y = 


which can be integrated once to yield 

,__J_ , A 

9 3 f f 2 ’ 


(32) 


with k a constant to be determined. A second equation and the determination of a comes 
from a leading order balance in the kinematic condition (10), which gives 


0 = 1, (s + /?£)/' + ( = (33) 

It can be seen that the leading order similarity solutions found here are identical to those 
predicted by the long wave model system. Solutions which are to be compared with the 
numerical solutions of the evolution equations are described next. 


4.3 Solution of the similarity equations. 

To fix matters we work wdth equations (25) and (26). Solutions must be obtained for — oo < 
£ < oo. The behavior for large £ is easily established from the equations to be 

m-K ih (33) 

This asymptotic behavior can also be deduced from the similarity transformations (23) since 
far away from the pinch region (i.e. as |£| -* oo) the solutions are expected to be independent 

of t ; it follows, then, that f(() ~ and g(Q ~ which make S 0 ~ \z\? and 

(Co ~ |r|~"A which are the outer solutions as \z\ — > 0. It follows from (34) that g vanishes as 
|£| tends to infinity which in turn implies that there is a point, £ 0 say, where g(( o) + = 0. 

Such a point is a removable singular point of equation (33) and requires a local analysis; £ 0 
can be shifted to the origin by use of the transformations 

f — > /, G(g) = g + j3 £ o y = i ~ £o- 

The point 7 = 0, then, is a removable singularity if 

(7(0) = 0, G"(0) = 2. 

This is a condition required by smoothness of solutions. A local analysis for (77 1 << 1 gives 
f{y) = /o + 7V2 + 7 4 .A + • ■ ■ ) G(rj) = 2i) + y 3 g 3 + . . . , 


where 


1 j = 3 + 2 0 

• /o_ 12(l +P)' 72(1 + fl) 2 ' 


(35) 
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with the remaining coefficients expressible in terms of a single parameter / 2 . 

To get a closed form solution defined implicitly it is useful to eliminate </(£) between (25) 
and (27), for example, to obtain 


r = iff (Lz UlMzl 

f (/+ 7^K/-/o)’ 


which can be integrated once to give 

L 


, 8+1 


(f _L 3+2/?.. V 

^ — m r + ■ a = An- 


k /(/ - M ' 12 


The substitution 


/ = 


1 


12(1 + j3) 

in (37) above, leads to the implicit solution 


— cosh 2 (0), 


( 12(1 + P )) 0 


f 

Jo 


cosh 6 


where 


A = 


f(n) 

12(1+ (3) 

f+k 


’ - dO 

— +Al). 

G(v) 

- LI 

2 + 13 


6(1 +13) 


cosh 2 («), 


(36) 


(37) 


(38) 


and ± corresponds to rj positive and negative respectively. Using the boundary conditions 
G'(-oo) = G'(Too) = dto we see that 


, _ i J2°(i/f)dn 

stfWPW 

which is seen to be identical to (28) once we note the symmetry of /. A numerical procedure 
is required to find admissible values of jS. The first, thing to note is that the ratio of the two 
integrals in (39) is independent of / 2 as can be seen by a rescaling of g in (38) (this has been 
established numerically also). Given a value of the scaling parameter f3, the solution /( tj ) is 
constructed from (38) by prescribing a value of 6 which implies an elevation / at an axial 
position 7/ determined by integrating (38) over the appropriate 0-range; with f(rj) known 
over an interval 0 < 77 < rj max , where T} max is sufficiently large for the asymptotic behavior 
(34) to be valid, an iteration in (3 is carried out to satisfy the eigenrelation (combining (35) 
and (39)): 



3 + 2/3 _ 1 JniAQffi, 

72(1 + /3) 2 3/ 0 ~(l/m' 
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These computations give a unique value — 0.175 correct to three decimals. Representative 
solutions are given in Figure 1 from which it is seen that the effect of is just a rescaling 
in 7] as expected. 

This Section has provided a fairly complete quantitative description of possible terminal 
pinching states of Stokes jets. The most notable feature of the analysis is that unique scal- 
ing exponents are fixed together with universal scaling functions to within a multiplicative 
constant which depends on the initial conditions. The uniqueness arises by solving an eigen- 
value problem for symmetric pinching profiles. The remainder of the article is concerned with 
the verification of the local analysis just described by detailed comparison with numerical 
solutions of the model system for different initial conditions. 


5 Numerical solutions. 


Numerical solutions have been obtained on axially periodic domains which can, without loss 
of generality, be normalized to have length 2 x. Since the pinching phenomenon is a local 
one the periodicity in the boundary conditions is not expected to play a fundamental role. 
The equations to be solved numerically are 

Sot T u'qSq, 4- —WozSo = 0, (41) 

(35oU?os 4 S 0 ) z = 0. (42) 


As will be shown shortly the initial condition 5o(0,s) = F(z) alone is required. It is clear 
from (41) that. 

(43) 

which provides a useful conserved quantity in controlling the accuracy of the computations. 
Equation (42) can be integrated once with respect to 2 to yield 


w Qz 


3 V So 2 So) ’ 


(44) 


where A(f) is a function of time to be determined. Integration of (44) with respect to z and 
use of periodicity yields a value for A(f) in terms of Sq: 


m = 


jS*(lfSo)dz . 

ti W {VSl)dz' 


(45) 


Equation (44) combined with (45) shows that an initial condition for S 0 alone is sufficient 
to determine the evolution. 

Spectral methods are used to solve the equations numerically, and all required integrations 
are also done spectrally. The result used is that the discrete fourier transform of a periodic 
function. q(z) say, and the integral of q(z) are related by the following expression which 
follows directly from the definition of fourier transforms: 



For example A (t) in (45) is computed by calculating the fourier transforms of (1/5 q) and 
(1/5 q) and forming the ratio of the first two fourier components. With these results available 
the numerical strategy is as follows: Equation (41 ) is marched fomard in time by specification 
of 5o(0, z ) alone; A is computed as indicated above at each level of a multi-level time scheme 
which then provides ivo z through (44) at that level; to update So in time using (41) the 
function iv o is also required. This is evaluated by transforming (44) into Fourier space by 
use of fast fourier transforms and inverting (ik)~ l H, where H is the discrete fourier transform 
of | — Jr-J. (note that the fourier transforms of (l/5o) and (1 / Sq) are already known 

at this stage from the computation of A). Finally the remaining derivative So z in (41) 
is evaluated pseudospect rally and the time integration is done in real space. A second 
order accurate predictor-corrector scheme as well as a third order Runge-Kutta method were 
implemented and tested. The computations described here were generated by the second 
order scheme. Certain symmetries in the equations were also utilized in the numerical 
work and to provide numerical solutions that can be directly compared with the similarity 
solutions described earlier. It can be seen that if So is an even function of z initially, then it 
will remain so for subsequent times since wq is then odd as seen from (44). More formally it 
is easy to show by integration by parts of (44) that 

< «'o >= ■ L (- < V >< *s 0 - 2 > + < V >< 2S 0 -' >) • (46) 

o < o o > 

where <(■)>= Jo [• )dz, and from which it follows immediately that w 0 is odd if So is even. 
With w o odd it is ensured that its zero fourier component vanishes and does not enter into the 
fourier inversion procedure of its computation outlined above. For general initial conditions, 
however, the inversion to find wo can only be performed if the first fourier mode (the one 
corresponding to k = 0) is provided. In the notation of (46) this is just < wo > which is 
easily found by performing the additional fourier transforms of z/Sq and z/Sq respectively, 
or quadrature methods. Both were implemented and used with no noticable change in the 
results. 

In our computations of symmetric solutions we use the following initial condition 

5(0, z) — a + bcos(z). (47) 

This choice gives a minimum in the initial condition at z = 7r (b > 0) and symmetry about 
this point with pinching first taking place there. The corresponding axial velocity is then 
zero at 2 = 7r and an odd function of 2 — 7r. Throughout the computation the evolution 
of < Sq > and < u'o > was monitored in order to confirm that the former is conserved 
and the latter is zero due to parity; in all the results given here these integral constraints 
varied only due to machine round-off errors. Besides the actual solutions at different times, 
it is important to monitor the evolution of A (t) and the value of the minimum as a function 
of t which we denote by S mtn (t). These quantities are crucial in our comparisons with the 
similarity solutions of Section 4. The axial position where the minimum is attained was also 
monitored in order to verify accuracy, since the minimum is stationary for the present choice 
of initial conditions. We note that the accuracy tests are crucial in determining acceptable 
solutions which can be used to verify the asymptotic theory especially since such comparisons 
are only meaningful when the jet radius is very small and the equations become singular. 
All solutions given here pass the integral accuracy tests outlined above. 
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Before presenting numerical results we consider the linear stability of the system (41) 
and (42) which is being solved numerically. Linearization is done about a constant value of 
S 0 taken to be a for consistency with the initial condition and about the zero state for iv 0 . 
Linearization of (41) and (42), then, and elimination of w 0zz by differentiation of the first 
equation with respect to z gives the following equation for the linear perturbation denoted 
by So also, 

Sozt — — 0. (48) 

DG 

The solution is readily found to be 

5 o (M) = So(0,r)e^. (49) 

The linear solution (49) is identical to the growth rate of the k = 0 mode of linear theory as 
expected (this follows from (1) and (2) once the non-dimensionalizations (3) are introduced). 
It can be seen from (49) that linear perturbations grow' in place, so that a depression (regions 
where ,?o(0, z) is negative) tends to grow reducing the local jet radius with the opposite 
happening in regions w r here there is a local elevation. Physically this implies that fluid 
is being pushed out of a necking region and towards bulging regions, a result w'hich is in 
line with the mechanism of capillary instability since surface tension tends to increase the 
pressure just below- the surface of a depression and decrease it near elevations, causing a 
fluid motion from high pressure regions to low pressure ones. i.From a numerical point of 
view the linear result (49) indicates that numerical short wave round-off error disturbances 
w'hich are inevitably introduced into the numerical simulation, are not subject to pathological 
instabilities as in Kelvin- Helmholtz or related problems (see Krasny [23] and Papageorgiou 
and Smith [24] for instance). 

The first set of numerical experiments has an initial condition with a — 0.5 and b = 0.1. 
The number of modes used is 512 and the time-step was 0.0005 by the end of the computation. 
The computation w r as terminated when the minimum jet radius, S m i n , became smaller than 
0.003. This happened at approximately t = 6.6. In order to achieve such a small value of S m i n 
the computational step was refined during the stage 6.5 < t < 6.6. Numerical convergence 
was checked by performing an identical computation with 256 modes; the only difference is 
that the higher resolution computation allows achievement of a slightly smaller minimum 
radius for a given accuracy. The evolution of the jet surface bo(t,z^ and the corresponding 
axial velocity iv 0 (t,z) are given in Figures 2(a), (b). Figure 2(a) depicts the evolution for 
0 < t < 6.5 while Figure 2(b) that during 6.5 < t < 6.6. As the Figures indicate the 
jet is pinching after a finite time; the radius vanishes first at z — ir. The final computed 
shape shown as a cross-section of the jet, along w'ith the evolution of the maximum value 
of S x over the spatial domain, is given in Figure 3. The jet. shape appears fairly flat near 
the pinch point (the slope is in fact zero there) and in particular the slope is bounded at 
all axial positions. This is an essential requirement for the validity of the long-wave model 
equations. The leading order axial velocity distribution has the following features: The flow' 
is stagnant (to leading order) at z = tt where a pinch first appears. Just to the right and 
left of the pinch point there are large axial velocities away from the pinch points; these axial 
velocities become infinite in magnitude as the singular time is approached and fluid drains 
out of the necking region at increasingly higher local speeds. We will use the results from 
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this numerical experiment to perform a detailed comparison with the asymptotic theory of 
Section 4, and in particular will derive the scaling exponents and the rate at which the force 
in the jet is going to zero as the singular time is approached. 

Figure 4 shows the evolution of S m i n (t), which is seen to approach zero after a finite 
time - the rate at which this is done is found to be linear (see later for an estimate of the 
slope), and as shown in the enlarged section the time of the singularity can be estimated by 
a least squares fit. The evolution of the maximum axial velocity with time, along with the 
corresponding evolution of the axial position where it is attained, is shown in Figure 5. These 
results strongly suggest that the axial velocity blows up at a finite time and this happens 
locally at some axial position. The rate at which the velocity blow's up will be compared 
with the self-similar solutions later. Note that due to symmetry there is a symmetrically 
placed minimum in it* o which blows up at the same rate as the maximum. Finally, in Figure 
6 we depict the evolution of A (t) given by (45) noting that physically this quantity represents 
the quasi-uniform force, to leading order, throughout the jet at different times. Clearly A (i) 
is approaching zero after a finite time; in fact for times larger than approximately 4.0 the 
dependence of A with t is established to be linear, a fact which is used to get an estimate of 
the singular time, t s , by a least squares fit. as above. 

The results just presented indicate a qualitative picture of the breakup phenomenon: the 
radius of the jet goes to zero linearly and so does the force in the jet. At the same time the 
axial velocity blows up after a finite time. These numerical results are used next to make a 
direct comparison with the similarity solutions constructed in Section 4. Since the pinching 
similarity solutions of the model equations are the same as those for the Stokes equations 
(see Sections 4.1 and 4.2 respectively) the comparison holds for both regimes. According to 
the ansatz (23), then, along with the solutions (38) it is easy to determine S mtn (t) given by 
the local similarity theory, 

Smin{t ) = {ts ~~ t) ^ ) (50) 

where ,6 = 0.175 (see Section 4.3). An expression for A (f) near the singular time follows 
similarly from the analysis of Section 4 and yields 

mo = m . -<> = 24( 1 + +% ) » ( 5l > 

where the expression (35) has been used to express k in terms of j3. The analytical ex- 
pressions (50) and (51) allow for a direct comparison between the theory and the numerical 
experiments. To carry this out, however, we require the value of the singular time, t s . This 
is not provided directly by the numerical solutions but can be estimated with tolerable ac- 
curacy from the data which make up Figures 4 and 6. A value of t s was obtained as follows: 
a least squares fit was applied to the data in Figure 4 for t > 5.0 and extrapolation was used 
to find t s the point where the curve (straight line) intersects the f-axis. This was repeated 

with least squares fits of smaller sets of data nearer the singular time with no change in the 

result. As a check the t s estimated by least squares fits of the \(t.) data (Figure 5) yields the 
same value of t s to four decimal places. This value is 

ts = 6.6405. (52) 
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Using this value of t a we can check the power law behavior of the computed solutions 
near the singular time. This is clone as follows: We take the sets of data (t, (U ■MOK 

(/, w max {t )) from Figures 4, 5 and 6 respectively for t > 5.9 (this is near enough the singularity 
for the asymptotic structures to provide a good approximation to the exact solutions) and 
produce plots of the sets {(log(f s — t), log^SGin^f))} and so on; the slopes of these lines 
should give the powers of r appearing in the similarity ansatz (23). The results are depicted 
in Figure 7; the slopes are estimated by least squares fits and are correct to the number of 
decimal places indicated on the diagrams. The results are also summarized in Table 1 below. 


Scaling law 

Asymptotic 

Numerical 

S m in{ts - 0 

1.0 

1.0 

A (t. ~ t) 

1.0 

1.0 

max{wo}(t s — t ) 

Cu 

1 

II 

1 

o 

GO 

CO 

-0.823 


Table 1 

The results just presented verify the power law behavior postulated by the similarity 
ansatz (23). The value of ,3 = 0.175 is therefore supported by the numerical solutions also, 
and next we make a further comparison of the multiplicative constants also. As shown 
above, the leading order behavior of S min {t) and A (t) according to the asymptotic theory of 
Section 4 as the singular time is approached is given by (50) and (51). In Figures 8(a), (b) 
we superimpose the numerically computed evolution of S min and A with the corresponding 
asymptotic forms (50) and (51). It is seen that at times larger than about 4.0 (which is still 
at least 2 time units from the singular time) agreement is excellent; the divergence of the 
two curves for smaller times is expected. 

The results of Figures 3-8 provide strong evidence that, the initial value problem of the 
model equations terminates in a singularity after a finite time according to the theoretical 
predictions set out in Section 4. The comparisons above were designed to verify the power 
law behavior of the solutions and in what follows we consider a comparison of the solutions 
of the model equation near the singular time with the scaling functions found from the 
asymptotic theory. The numerical experiment we have been concentrating on terminates at 
a singular time estimated to be t s = 6.6405. We describe next how to construct the scaling 
functions /(£) and g(() from our numerical solutions near the singular time: Given data 
S 0 (M) and w 0 {t,z) it is easy to compute £ and the corresponding values of /(£) and g(() 
by direct substitution of r = (i s — f) into the forms (23). Different scaling functions are 
obtained for different r but all results should collapse to universal curves as r — > 0. This is 
indeed the case as depicted in Figures 9(a), (b) corresponding to /(£) and #(£) respectively. 

5.1 General initial conditions; non-parity solutions 

The results described so far have been computed by imposing a symmetry on the evolution 
equations, restricting solutions to shapes which are symmetric and axial velocities which are 
anti-symmetric about z — tt. This was done in order to obtain an accurate set of data to use 
in the evaluation of the asymptotic self-similar theory. The self-similar theory, however, is a 
local one and it does not require the parity requirements mentioned above for all time. We 
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would expect locally symmetric solutions, then, to appear as correct leading order dynamics 
in the breakup described by the model equations (41) and (42) starting from general initial 
conditions which do not ensure parity for subsequent times. Equations (41) and (42) were 
solved numerically (see earlier) using n = 512 modes and the initial condition 

r<o(x) = 0.5 + 0.1(sin(.r) + cos(.t)). 

The evolution up to t = 5.0 of the interfacial shape and the corresponding axial velocity 
is given in Figure 10. Due to the loss of symmetry, the jet first breaks at a point not 
equal to tt. The results of this numerical experiment are used next to confirm the local 
validity of the self-similar solutions. In Figure 11 we give the evolution of 5 m ,„(f) and A(f) 
up to the time when S mtn is less than 0.005. Superimposed on these curves are the results 
given by the asymptotic forms (50) and (51) with ,3 determined by the similarity solutions 
({3 - 0.175) and t s determined numerically from the data, by first confirming that the curves 
are straight lines as t —> t s , followed by extrapolation to obtain t s . It is seen that agreement 
is excellent. A more severe test is a direct numerical check of local symmetry. This is done 
as follows: The profile (jet shape) is taken from the last computed time station and the 
minimum point on the curve is located; this gives the axial position where the jet radius is 
at its smallest; if we denote this position by .ro, our objetctive is to show that the jet shape 
is symmetric (at least locally) about Xq. A good way to see this graphically is to plot the 
jet shape for x > .ro and superimposed onto this to plot the shape computed for x < .To 
but reflected about the line x — Xo. If there is local symmetry about x — x 0 then the two 
curves will coincide for a range of axial positions in the neighborhood of To- The results are 
shown in Figure 12(a) with the circles denoting the reflected shape. We see that the two 
curves are indistinguishable for a large range of values about x 0 and so we conclude that 
the solution at pinching is locally symmetric. The corresponding axial velocity distribution 
becomes asymetric as a pinch is formed, and the confirmation of this based on our numerical 
results is given in Figure 12(b) which depicts plots of (x,Wo(x)) for x > .r 0 plotted along 
with (.t, —wo(xo-x)). Again the curves are indistinguishable confirming the local self-similar 
theory. Our numerical solutions, therefore, provide strong evidence that the local self-similar 
solutions described in Section 4 are robust in that they provide a local theory with symmetry 
independent of initial conditions. The initial conditions affect two things, however: (i) The 
singular time, and, (ii) the scaling function (and thus the final shape) at breakup. These 
two features are studied in more detail next by carrying out a parameter study as the initial 
condition varies. 


6 Breakup times and scaling functions for different 
initial conditions. 

It is expected that different initial conditions provide different breakup times. For example, 
it would be expected that the breakup time decreases as the minimum initial jet radius is 
decreased. In what follows we try to quantify such statements by carrying out a series of 
numerical experiments - this is done for both symmetric and asymmetric solutions. Before 
presenting results we define what we mean by breakup time. All runs (unless otherwise 
stated) are followed up to times when the minimum jet radius first becomes less than 0.005. 
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The evolution of A (t) is then used to obtain an estimate of the singular time by extrapolation 
after utilizing the linear form of A (t) near the singular time (see Section 5). The self-similar 
shapes at breakup are then approximated by the data from the last computed time, after 
applying the transformations of (23). 

One aim of carrying out such extensive computations is in the evaluation of linear stability 
theory, estimates of the breakup time for instance, as compared to the nonlinear dynamics. 
Working within the framework of the one-dimensional model, we can use linear theory to 
obtain an empirical estimate of the breakup time. As shown in (49) linear solutions grow 
exponentially at a rate exp(f/6a); the total solution which results from an initial condition 
h{z) is 


S(z , i) = a + 8h(z) exp(t/6a), 


(53) 


where 8 is the infinitessimally small amplitude of the perturbation as is usual in linear 
theories. Denoting the minimum of h(z) over the domain by —h 0 < 0, we can use (53) to 
predict a time, t L say, when S(zJl) first becomes zero. This is easily calculated to be 

(l = 6<,ln (^)’ (54) 

and represents the breakup time as predicted by linear theory. We emphasize that this 
calculation is empirical in that, for example, t L is a time outside the range of validity of 
linear theory (linear theory is valid when S(z,t) is near r = a). Similar ideas have been used 
in [25] in describing the rupture of free viscous films in the presence of van der Waals forces. 


6.1 Symmetric initial conditions 

Symmetry in the jet shape is preserved throughout the evolution if the initial condition is 
s y jyj metric. In the results that follow we computed breakup times as a function of an initial 
amplitude ej where 

S'o(2,0) = 0.5 + Ci cos(.r). 

The parameter £i is also a measure of the initial energy provided by the disturbance. We 
note that e, < 0.5, otherwise the initial condition has zero radius to start with - the breakup 
time in this case is defined to be zero. The results are summarized in Figure 13 which 
depicts the variation of the breakup time, t s , with t\. The figure also contains the breakup 
time t L predicted by the linear result (54), noting that in this case h 0 = 1. It is seen 
that agreement between the two is only reasonable when t\ is small as would be expected 
(we emphasize once more, however, that the linear estimate albeit empirical can be useful 
in some instances). For example, an initial perturbation ej = 0.05 (a 10% disturbance of 
the unperturbed jet radius) gives a breakup time t s = 9.04 while linear theory predicts 
t L = 6.91, an underestimate of approximately 26%. At a 20% initial perturbation the error 
is approximately 27%. An additional feature of the numerical results, which is in line with 
intuition, is that the breakup time decreases as the initial amplitude increases. In fact for 
initial amplitudes larger than about 0.2 (i.e. a 40% perturbation of the unperturbed radius) 
the breakup time varies almost linearly with ex- The numerical values that make up Figure 
13 are also given in Table 2 below. 
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Cl 

.001 

.005 

.01 

.02 

.03 

.04 

.05 

.075 

.1 

.125 

.15 

.175 

is 

21.13 

16.28 

14.16 

12.00 

10.71 

9.78 

9.04 

7.66 

6.64 

5.82 

5.13 

4.53 

tl 

18.64 

13.82 

11.74 

9.66 

8.44 

7.58 

6.91 

5.69 

4.83 

4.16 

3.61 

3.15 

Cl 

.2 

.225 

.25 

.275 

.3 

.325 

.35 

.375 

.4 

.425 

.45 

.475 

.5 

ts 

4.00 

3.52 

3.08 


2.31 

1.95 

1.63 

1.32 

1.03 

.751 

.487 

.230 

0 

tL 

2.75 

2.40 

2.08 

1.79 

1.53 

1.29 

1.07 

0.86 

0.67 

0.49 

0.32 

0.15 

0.0 


Table 2 

6.2 Non-symmetric initial conditions 

The results that follow were generated from an initial condition of the form 

5 o (z,0) = 0.5 + ei(sin(x) + cos(.t)). 

Noting that this can be re-written as 

S o {z,0) = 0.5 + cos(.t — tt/4), 

we see that the value e\ = \/2/4 ss 0.35355 is an upper bound for the initial perturbation 
amplitude. The breakup times were computed as before, and the estimate provided by 
equation (54) is valid with 6 = \/2ei now. Figure 14 shows the variation of the computed 
breakup time with the initial amplitude level ei, together with the linear result ti. Table 3 
below provides the numerical values that comprise the figure. 


Cl 

.001 

.005 

.01 

.02 

.03 

.04 

.05 

.075 

.1 

ts 

20.09 

15.22 

13.08 

10.90 

9.59 

8.63 

7.86 

6.43 

5.36 

tL 

17.60 

12.78 

10.70 

8.62 

7.40 

6.54 

5.87 

4.65 

3.79 

Cl 

.15 

.2 

.25 

.275 

.3 

.31 

.32 

.33 

.34 

.35 

ts 

3.76 

2.56 

1.58 

1.16 

.76 

.61 

.46 

.32 

.18 

.05 

tL 

2.57 

1.71 

1.04 

0.75 

0.49 

0.39 

0.30 

0.21 

0.12 

0.03 


Table 3 


Qualitatively, then, the behavior is similar to the symmetric case. A more direct compar- 
ison between the two types of numerical experiment is made in Figure 15. Here the breakup 
times are plotted as a function of the energy of the initial perurbation defined by 




Applying this to the initial conditions used above gives the following expressions for the 
perturbation energy corresponding to symmetric, E s say, and asymetric, E a say, initial 
conditions, 



Figure 15 shows collectively the breakup times as a function of E s and E a ■ It can be 
concluded from these results that the jet breaks sooner for the case of asymmetric initial 


17 





conditions; for the particular choice of initial conditions used here the breakup times for the 
symmetric case are always larger than the corresponding ones for asymmetric conditions, at 
a given equal initial perturbation energy. Asymmetry appears to enhance breakup, then. 

6.3 Scaling functions for different initial conditions 

Here we construct the behavior of the solutions near the singular time for different initial 
conditions. For brevity we consider pinching solutions computed from symmetric initial in- 
terfacial elevations - non-symmetric conditions are treated in the same way and have already 
been shown to provide locally symmetric solutions at breakup (see Figure 12). The main 
objective of this section is to provide numerical evidence that different initial conditions will 
pinch according to a single scaling function when suitably normalized. We constructed scal- 
ing functions near pinching by solving the initial value problem for different Ci in (6.1). The 
scaling functions were constructed as described earlier (for instance the methods used in the 
construction of Figure 9) and by use of the computed estimates of singular times from Table 
2. In what, follows we concentrate on the interfacial shape scaling functions derived from 
numerical solutions starting from initial conditions with ei = 0.005,0.1,0.2,0.3. According 
to the similarity solutions of Sect ion 4, all scaling functions coincide at £ = 0 where they take 
the value 1/(12(1 + /?)). This will be exhibited in the construction of the scaling functions 
below, but can also be seen by consideration of the variation of u min (t ) as the singular time 
is approached - the variation of u min (f) with t near the singular time t s for a given initial con- 
dition, should be linear with slope 1/(12(1 + 0)) (see (50)). Figure 16 provides the variation 
of u mtn (t) for different initial conditions d = 0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3 labeled 
on the figure and superimposed with these numerical results is the asymptotic behavior near 
the singular time (see above and equation (50)) shown wuth open circles. It is seen that 
agreement is achieved near the singular times as expected. 

The scaling functions for different initial conditions are considered next. As showm by 
the results of Figure 16, all scaling functions are equal to the universal constant 1/ (12(1 + 
/3) ^ 0.0709 at £ = 0. Numerically constructed scaling functions from initial conditions 
characterized by the amplitudes fi = 0.005,0.1,0.2,0.3 are shown in Figure 17(a) which 
indicates the self-similar nature of the different terminal scaling functions. It can be seen from 
the analytical solution (38) that the only difference between scaling functions is expected to 
appear through the constant A in (38) which in turn depends on initial conditions; the role 
of this constant is to stretch the axial coordinate ?/ by different amounts for different initial 
conditions. Our numerical simulations and singular states fully support this property of the 
solutions as is demonstrated next: Choose the computed scaling function corresponding to 
fl = 0.1 as the reference function, foijj) sav. According to the theory, for each of the other 
initial conditions shown in Figure 17(a), a number, c say, can be found so that the change 
of variables rj — ► c?/ maps the given scaling function into the normalizing one. The number 
c is different for different initial conditions and was computed by calculating the ratio T] 0 /r), 
where fo{rjo) = hiVi) = h with h a fixed interfacial amplitude - different values of h produce 
the same ratios as expected (the value of h had to be found by interpolation due to the non- 
uniformity of the grids). This procedure, then, enables all scaling functions to be collapsed 
onto fo(i]). The result of this calculation is given in Figure 17(b) which provides additional 
evidence for the validity of the asymptotic theory near the singular time. 
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7 Computation of satellite drops 

In previous Sections we provided numerical and analytical evidence of the form of pinching 
solutions admitted by the model evolution equations. Even though the equations are rela- 
tively simple compared with the full governing system (the Stokes equations plus nonlinear 
interfacial conditions) it has been established that the model system captures a lot of the 
nonlinear stages of the evolution leading to pinching. In this Section we use the model to 
compute final breakup states which are remeniscent of the observed phenomenon of mother 
and satellite drop formation. To compute such solutions the initial conditions need to be 
chosen appropriately. In order to illustrate things we concentrate on symmetric solutions, 
and in particular we use a one-parameter family of initial conditions given by 

S o (2,0) = 0.3 + 0.05cos(z) + e 2 exp(-(7r - zf). (55) 

This initial condition is symmetric, periodic in z and has a localized hump at z — tt over the 
longer wavelength depression provided by the cos(^r) perturbation. The advantage of such 
an initial condition is that the jet pinches after a finite time with the radius vanishing simul- 
taneously at tw 7 o distinct axial positions. Just beyound pinching, then, a system of mother 
and satellite drops emerges. Sample results are depicted in Figure 18 for five successively 
increasing values of e 2 . The Figure shows the cross section of the jet just before pinching; 
the computational domain is 27r-periodic in the axial direction and for better visualization 
the solution has been dipicted to include four periods. We have used the computational 
methods described in earlier methods to verify that the pinching takes place according to 
the self-similar solutions given in Section 4, and we have computed singular times as w 7 ell as 
the volume of the satellite drops. These are estimated by evaluating the integral 



where and z+ are the minima in So to the left and right of the satellite drop respectively. 
The integral lo = x S * S$dz is a conserved quantity which is equal to the total volume of 
the jet at / = 0. The ratio V s /Vq provides a measure of the size of the drops which form 
after pinching. The present numerical experiments yield the results given below r in Table 4, 
and summarized in Figure 19. 
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Table 4 

It can be seen from the results that relatively small changes in the initial conditions can 
lead to large changes in the final drop size. Note that the initial volume of the jet increases 
by approximately 3% when e 2 is increased from 0.03 to 0.05, w'hile the volume of the satellite 
drop formed increases by more than 100%. It appears, then, that- small increases in e 2 enable 
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the satellite drop to drain fluid from the adjoining main jet yielding a wide range of drop 
sizes. At the same time the jet breaks sooner as seen from the variation of t s . The trends 
depicted in Figure 19 will be studied in more detail by more extensive computational searches 
in the phase space of initial conditions. 


8 Conclusions 

An asymptotic theory has been used to derive a system of nonlinear evolution equations to 
model the dynamics of viscous fluid jets under the action of capillary forces. The theory 
is a long wave expansion which allows for fully nonlinear interfacial amplitudes and the 
equations can predict pinching. Theoretical descriptions of the pinching have been given by 
use of similarity solutions valid as the jet radius tends to zero after a finite time. The scaling 
exponents are unique and the scaling functions form a one parameter family of similarity 
solutions depending on the initial conditions. The predictions of the asymptotic solutions 
near pinching have been confirmed, with excellent agreement, by extensive direct simulations 
of the initial boundary value problem. The solutions found here are possible terminal local 
states of the full Stokes equations and can be used to benchmark the accuracy of direct 
simulations (e.g. boundary integral or boundary element techniques), as well as to provide 
correct initial conditions to continue the computations beyond pinching. Such analyses have 
been described elsewhere for viscous jets (see Papageorgiou [26]). These comments are also 
relevant to viscous jets possessing inertia and modelled by an analogous system of evolution 
equations (see Appendix A). 

In addition, it has been shown numerically that the model can produce pinching solutions 
with the radius of the jet vanishing at the same time at two distinct axial locations; this 
heralds the formation of a mother-satellite drop system and the model has been used to gain 
quantitative information on the distribution of satellite drop volumes as a function of initial 
conditions (or equivalently initial perturbation energy). A considerable amount of numerical 
experiments remain to be carried out in order to obtain overall trends in both the Stokes 
and Navier-Stokes models. 
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Appendices 


A Navier-Stokes regimes 

Using the capillary scales (3) for dimensionless quantities and starting from the Navier-Stokes 
equations one is led to the following system: 


u 

R t {u t + uu r + u IV 2 ) = —p r + Art 

(56) 

R e (w t + uw r + ww.) = —p z + Arc, 

(57) 

~(ru) r + r o z = 0, 
r 

(58) 


where the Reynolds number R t = pp- (1 / R e is the familiar capillary number). The boundary 
interfacial conditions are the same as (8)-(10). Proceeding as in Section 3 we introduce a long 
axial length-sacle by the transformation -p — » ep. Our objective is to find an asymptotic 
solution in powers of e 2 which includes, to leading order, the effects of unsteadiness and 
nonlinearity in the bulk motion. In particular we wish to increase R e from its zero value for 
Stokes flows, to a value when unsteadiness and nonlinearity first enter to yield a canonical 
evolution system. The appropriate expansions for the flow parameters are those of Section 
3, see equations (11 )-(14), together with the scale 

R e = c 2 k, k = 0(1). 

The axial momentum equation (57) gives to the first two orders 

wo = w 0 {z,t), (59) 

k{u'o t + U’oU’Or) = ~POz + lt ? lrr H W’lr + WQz 2 -s (60) 

while the radial momentum equation and the continuity equation give, to leading order, 

~~P0r + UOrr H UOr T = 0, (61) 

r r l 

u 0 = --riv 0 z- (62) 

It follows by substitution of (62) into (61) that 

POr = 0 => Po=Po(z,t). 

This enables the calculation of po from leading order terms in the normal stress balance 
condition (9), and yields 

Po = Jr - i»o*. (63) 

‘-’0 

We are now r in a position to derive the first evolution equation from the tangential stress 
balance equation (8). This is done as follows: The velocity correction is easily evaluated 
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from (60) since the forcings w 0 and p 0 are independent of r. This solution, along with the ex- 
pression (62) for M 0 are substituted into (8) to yield the following leading order contributions 
after evaluation on r = So(z,t ): 


U’Or = o, 

— ^5oit’o; + vho [k(u’o< + n’o ?, ’o c ) — 2wq zz -f (1 /*S'o)z] — 3.S'o,t(,’oz = 0. 
Z Z 


(64) 

(65) 


A second evolution equation comes from the leading order contribution of the kinematic 
condition (10). Re-arrangement of (65) together with the kinematic contribution gives the 


, 3(5 0 2 u , o,) i 

K\Wot + W 0 W 0z ) = q 2 


1 

2‘ 


and So, 


( 1 ^ 

1 c j ’ 

wo / Z 

(66) 

& 

hi 

II 

p 

(67) 


The order one parameter k can be scaled out of the problem by a change of variables and 
so we can take k = 1. 

Equations (66) and (67) have been derived by Eggers and Dupont [17] by using a Tay- 
lor expansion in the r coordinate, and by the present approach by Papageorgiou [26] in a 
slightly different regime. The present approach is a systematic single parameter asymptotic 
expansion with all terms appearing in the equations being of the same order in e. Equations 
(66) and (67) admit, pinching similarity solutions (see [18], [19] and [26]) of the form 

So — (is ~ u ’o — (i* ~ i) t 9(0i £ ~ (j s _ /)i/2’ 


where t s > 0 is the time of pinching. The scaling functions appearing in (68) satisfy 

3/ + 41 + 6^21 = 1(3 + &')+!»', (69) 

(3 + M2)/'-(l-3'/2)/ = 0, (70) 

with primes denoting ^-derivatives. Eggers [18] and [19] has given numerical solutions of 
these scaling functions which do not depend on arbitrary constants and so describe pinching 
states which are independent of initial conditions. W'e note also that the system (69) and 
(70) arises by looking directly for pinching solutions of the Navier-Stokes equations without 
first deriving the intermediate asymptotic equations by introduction of e. Instead, the small 
parameter which allows for an asymptotic development is the smallness of t , t near the time 
of pinching. This analysis (as well as the equivalent one for Stokes flows) has been carried 
out by Papageorgiou [15] where the dynamics beyond pinching where also modelled by use of 
overall mass and moment um balances. Further numerical studies of (66) and (67) are under 
way to evaluate the model for such phenomena as satellite drop formation and comparisons 
with empirical breakup times obtained from linear theory, as well as direct comparisons with 

experiments. 
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List of Figures 


Figure 1 Solutions of the similarity equations for different values of f? shown on the 
figure; left - interfacial shape, right - axial velocity. 

Figure 2 Solution of the evolution equations for symmetric initial conditions: (a) 0 < 
t < 6.5, (b) 6.5 < t < 6.6. The computations were stopped when S mtn < 0.003. 

Figure 3 (a) The final computed jet shape at t = 6.6 when 5 mm = .0028; the figure 
represents a cross-section of the jet in a plane containing the jet axis, (b) Evolution of 
the largest magnitude of the interfacial slope - boundedness of this is crucial for the long 
wave/slender jet theory to be valid. 

Figure 4 Evolution of the minimum jet radius for symmetric initial conditions. The 
enlargement shows the linearity as t — » t s and the extrapolation to compute an estimate for 
t s . 

Figure 5 Symmetric initial conditions, (a) Evolution of the maximum axial velocity, (b) 
Evolution of the axial position where the axial velocity attains its maximum. 

Figure 6 Symmetric initial conditions. Evolution of A (/) (cf equation (45)), the leading 
order component of the force in the jet. The enlargement show's the linearity as t — ► t s . 

Figure 7 Symmetric initial conditions. Log-log plots estimating the pow r er law behavior 
of (a) * (b) A (t), and, (c) U’omaxfO? f > ^s- 

Figure 8 Symmetric initial conditions. Comparison of the similarity solutions with the 
computed ones, (a) S mt „(<), (b) A (t). The curves labeled asymptotic correspond to the 
theoretical predictions. 

Figure 9 Symmetric initial conditions. Convergence to scaling functions, (a) Interfacial 
shape in the pinch region, (b) Axial velocity in the pinch region. 

Figure 10 Solution of the evolution equations for asymmetric initial conditions, (a) 
Evolution of the jet shape, (b) evolution of the axial velocity. 

Figure 11 Asymmetric initial conditions. Evolution of S m i n (t), (a), and A (f), (b), to- 
gether with theoretical predictions. 

Figure 12 Asymmetric initial conditions. Verification of the symmetry of the interface, 
(a), and the asymmetry of the axial velocity, (b), near the pinch point. The solid lines 
denote the solution to the right of the pinch point, and the circles denote that to the left, 
after reflection about the pinch point. 

Figure 13 Symmetric initial conditions. Breakup times for different initial amplitudes, 
o - numerical calculation, * - prediction of linear theory. 

Figure 14 Asymmetric initial conditions. Breakup times for different initial amplitudes, 
o - numerical calculation, * - prediction of linear theory. 

Figure 15 Breakup times for symmetric and asymmetric initial conditions plotted as 
functions of the initial perturbation energy, o - symmetric, * - asymmetric. 

Figure 16 Variation of 5 mtn (0 for different symmetric initial conditions, ei = 0.005, 
0.01, 0.02, 0.03, 0.05, 0.1, 0.2, 0.3. The open circles show the asymptotic result according to 
equation (50) and the magnitude of their slopes is equal to 1/(12(1 + /?)). 

Figure 17 Symmetric initial conditions. The interfacial scaling functions near pinching 
for different initial conditions, (a) Non-normalized, (b) normalized with the scaling function 
corresponding to ei = 0.1. 
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Figure 18 Runs leading to satellatite drop formation. Numerical solutions generated for 
different initial conditions (55) and e 2 = 0.03,0.035,0.04,0.045,0.05 shown on the Figure. A 
jet cross section is shown just before breakup. 

Figure 19 Satellite drop formation, (a) Variation of the ratio of satellite drop volume 
to total jet volume with initial conditions; (b) Corresponding jet breakup times. 
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General data; comparison of similarity solutions with computation 
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Breakup times: Symmetric initial conditions 



Figure 13 
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Breakup times: Non-symmetric initial conditions 



Figure 14 
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Comparison of breakup times for symmetric and asymmetric 1C 
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Non-normalized interfacial scaling functions for different initial conditions 
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Normalized interfacial scaling functions for different initial conditions 
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